home *** CD-ROM | disk | FTP | other *** search
/ Die Ultimative Software-P…i Collection 1996 & 1997 / Die Ultimative Software-Pakete CD-ROM fur Atari Collection 1996 & 1997.iso / g / gnu_c / gpplib22.zoo / libsrc / xnormal.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1993-11-30  |  1.7 KB  |  61 lines
  1. /* 
  2. Copyright (C) 1988 Free Software Foundation
  3.     written by Dirk Grunwald (grunwald@cs.uiuc.edu)
  4.  
  5. This file is part of the GNU C++ Library.  This library is free
  6. software; you can redistribute it and/or modify it under the terms of
  7. the GNU Library General Public License as published by the Free
  8. Software Foundation; either version 2 of the License, or (at your
  9. option) any later version.  This library is distributed in the hope
  10. that it will be useful, but WITHOUT ANY WARRANTY; without even the
  11. implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  12. PURPOSE.  See the GNU Library General Public License for more details.
  13. You should have received a copy of the GNU Library General Public
  14. License along with this library; if not, write to the Free Software
  15. Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  16. */
  17. #ifdef __GNUG__
  18. #pragma implementation
  19. #endif
  20. #include <builtin.h>
  21. #include <xrandom.h>
  22. #include <xnormal.h>
  23. //
  24. //    See Simulation, Modelling & Analysis by Law & Kelton, pp259
  25. //
  26. //    This is the ``polar'' method.
  27. // 
  28.  
  29. double Normal::operator()()
  30. {
  31.     
  32.     if (haveCachedNormal == 1) {
  33.     haveCachedNormal = 0;
  34.     return(cachedNormal * pStdDev + pMean );
  35.     } else {
  36.     
  37.     for(;;) {
  38.         double u1 = pGenerator -> asDouble();
  39.         double u2 = pGenerator -> asDouble();
  40.         double v1 = 2 * u1 - 1;
  41.         double v2 = 2 * u2 - 1;
  42.         double w = (v1 * v1) + (v2 * v2);
  43.         
  44. //
  45. //    We actually generate two IID normal distribution variables.
  46. //    We cache the one & return the other.
  47. // 
  48.         if (w <= 1) {
  49.         double y = sqrt( (-2 * log(w)) / w);
  50.         double x1 = v1 * y;
  51.         double x2 = v2 * y;
  52.         
  53.         haveCachedNormal = 1;
  54.         cachedNormal = x2;
  55.         return(x1 * pStdDev + pMean);
  56.         }
  57.     }
  58.     }
  59. }
  60.  
  61.